/*************************************************************************
 * 
 *  fill histogramms at the end of an event
 * 
 *************************************************************************/

#include "EventAction.hh"

using namespace std;
using namespace CLHEP;

EventAction::EventAction(HistoManager *HM, Data* Data)
  : MiniBallEventAction(HM), histoManager(HM), data(Data)
{
  
  collectionID.resize(data->NumberOfHC);
  for(int i=0;i<data->NumberOfHC;i++)
    collectionID[i] = -1;

  RingNumber.resize(data->settings->NumberOfCDSectors);
  RingEnergy.resize(data->settings->NumberOfCDSectors);
  StripNumber.resize(data->settings->NumberOfCDSectors);
  StripEnergy.resize(data->settings->NumberOfCDSectors);
  Edet.resize(data->settings->NumberOfCDSectors);
  
}

EventAction::~EventAction()
{
}

void EventAction::BeginOfEventAction(const G4Event*)
{
  if(data->VerboseLevel > 4)
    cout<<"EventAction::BeginOfEventAction start"<<endl;

  for(int i=0;i<data->NumberOfHC;i++)
    {
      if(collectionID[i]<0)
	{
	  if(data->VerboseLevel > 4)
	    cout<<"trying to find "<<i<<". HitCollection "<<data->HCName[i]<<" ... "<<flush;
	  collectionID[i] = G4SDManager::GetSDMpointer()->GetCollectionID(data->HCName[i]);
	  if(data->VerboseLevel > 4)
	    cout<<collectionID[i]<<endl;
	}
    }

  if(data->VerboseLevel > 4)
    cout<<"EventAction::BeginOfEventAction stop"<<endl;
}

void EventAction::EndOfEventAction(const G4Event* evt)
{
  if(data->VerboseLevel > 5)
    cout<<"EventAction::EndOfEventAction start"<<endl;
  
  //cout << "data->Reaction " << data->Reaction << endl;//kwimmer
  bool gammafound = false;

  if(histoManager){
    
    G4int event_id = evt->GetEventID();
    G4HCofThisEvent * HCE = evt->GetHCofThisEvent();

    size_t nb_detectors = histoManager->GetNbOfDetectors();

    for(size_t det_id = 0; det_id < nb_detectors; det_id++){
      //every detector has its own HitsCollection 
      G4int HCID = histoManager->GetMBHitsCollectionID(det_id);
		if(HCID < 0){
	G4cerr << "MiniBallHitsCollection of detector id "<< det_id << " is missing" << G4endl;
	continue;
      }
      MiniBallHitsCollection* HPGeHC = (MiniBallHitsCollection*)(HCE->GetHC(HCID));
      if(HPGeHC != NULL){
	//Number of hits in the current detector during this event
	G4int n_hit = HPGeHC->entries();
	
	for(G4int hit_nb=0; hit_nb<n_hit; hit_nb++){
	  //every hit is prozessed singularly
	  histoManager->StartOfInteraction();
	  
	  //Energie deposited in the detector
	  G4double Edep = (*HPGeHC)[hit_nb]->GetEdep();
	  if(Edep > 0){
	    gammafound=true;
	    //position in detector
	    G4ThreeVector positionID = (*HPGeHC)[hit_nb]->GetPosInDet();
	    //position in world
	    G4ThreeVector positionIW = (*HPGeHC)[hit_nb]->GetPos();
	    G4int trace_id = (*HPGeHC)[hit_nb]->GetTraceID();
	    
	    size_t det_id = (*HPGeHC)[hit_nb]->GetDetectorID();
	    size_t seg_id = (*HPGeHC)[hit_nb]->GetSegmentID();


	    histoManager->SetMBInteraction(det_id, seg_id, Edep, event_id, trace_id, positionID, positionIW);
	  }
	}
      }//if(HPGeHC)
      
    }//for det_id
  }//if histomanager
  fill = false;
  if(data->Coulex&&data->DownScale>0){
    double tmp = RandFlat::shoot(0.,(double)data->DownScale); 
    if(data->Reaction==0||gammafound){
      fill = true;
    }
    else if(tmp<1){
      fill = true;
    }
  }
  else{
    fill = true;
  }


  for(i = 0;i < data->NumberOfHC;i++){
	 if(data->VerboseLevel > 5)
      cout<<"Starting on HitCollection nr. "<<i<<endl;
	 if(histoManager && collectionID[i] > -1 && !data->Efficiency){
	  
      G4HCofThisEvent* HCE = evt->GetHCofThisEvent();
      HitsCollection* FirstHC = NULL;
      HitsCollection* SecondHC = NULL;
      if(HCE&&fill){
	if(data->VerboseLevel > 5)
	  cout<<"Working on HitCollection "<<data->HCName[i]<<" (nr. "<<i<<")"<<endl;

	//cout << "hit in detector type " << data->DetectorType(i) << endl;
	switch(data->DetectorType(i)){
	case 0: // forward barrel	  
	  FirstHC = (HitsCollection*)(HCE->GetHC(collectionID[i])); //FdeltaE
	  SecondHC = (HitsCollection*)(HCE->GetHC(collectionID[i+1])); // FE
	  BD(true, FirstHC, SecondHC);
	  i++;		      
	  break;
	case 1: // backward barrel
	  FirstHC = (HitsCollection*)(HCE->GetHC(collectionID[i])); //BdeltaE
	  if(data->BarrelType==0)
	    BD(false, FirstHC);
	  else{
	    SecondHC = (HitsCollection*)(HCE->GetHC(collectionID[i+1])); // BE
	    BD(false, FirstHC, SecondHC);
	    i++;		
	  }
	  break;
	case 2: //FCD
	  FirstHC = (HitsCollection*)(HCE->GetHC(collectionID[i]));
	  if(data->CDPad){
		 SecondHC = (HitsCollection*)(HCE->GetHC(collectionID[i+1]));
	    CD(true, FirstHC, SecondHC);
	    i++;
	  }
	  else
	    CD(true, FirstHC);
	  break;
	case 3: //BCD
	  FirstHC = (HitsCollection*)(HCE->GetHC(collectionID[i]));
	  if(data->CDPad){
	    SecondHC = (HitsCollection*)(HCE->GetHC(collectionID[i+1]));
	    CD(false, FirstHC, SecondHC);
	    i++;
	  }
	  else
	    CD(false, FirstHC);
	  break;
	default:
	  break;
	}//switch(data->DetectorType(i))
      }//if(HCE)
    }//if(histoManager && collectionID[i] > -1 && !data->Efficiency)
  }//loop over HCs  

  histoManager->EndOfEvent();
  if(data->VerboseLevel > 5)
    cout<<"EventAction::EndOfEventAction stop"<<endl;
}

void EventAction::CD(bool forward, HitsCollection* FirstHC, HitsCollection* SecondHC){
  if(data->VerboseLevel > 2)
    cout<<"Starting CD Hits collection, clearing vectors"<<endl;

  int quads = data->settings->NumberOfCDSectors;
  for(int j=0;j<quads;j++){
    RingNumber.at(j).clear();
    RingEnergy.at(j).clear();
    StripNumber.at(j).clear();
    StripEnergy.at(j).clear();
    Edet.at(j) =0;
  }
  if(data->VerboseLevel > 2)
    cout<<"clear done"<<endl;

  int k;
  double energy;
  vector<double> Etot(quads);

  for(int j=0;j<FirstHC->entries();j++){
    energy = (*FirstHC)[j]->GetEdep()/histoManager->GetEnergyUnit();
    
	 int quadr = (*FirstHC)[j]->GetSiSegmentID();
    if(quadr >= 0 || data->WriteAll){ 
      //SegmentID is only set for hits in the active area of the CD
      if(quadr >= quads){
	cerr<<"error in hitcollection "<<i<<" hit "<<j<<" segmentid = "<<quadr<< " > " << quads <<endl;
	continue;
      }
      Etot[quadr]+=energy;
      //check whether this ring was already hit
      for(k=0; k<(int)RingNumber[quadr].size(); k++){
	if(RingNumber[quadr][k] == (*FirstHC)[j]->GetAnnularID()){
	  RingEnergy[quadr][k] += energy;
	  break;
	}
      }
      if(k == (int) RingNumber[quadr].size()){
	RingNumber[quadr].push_back((*FirstHC)[j]->GetAnnularID());
	RingEnergy[quadr].push_back(energy);
	if(data->VerboseLevel > 5){
	  cout<<"adding hit to new ring!"<<endl;
	}
      }
      //check whether this strip was already hit
      for(k=0; k<(int)StripNumber[quadr].size(); k++){
	if(StripNumber[quadr][k] == (*FirstHC)[j]->GetRadialID()){
	  StripEnergy[quadr][k] += energy;
	  break;
	}
      }
      if(k == (int) StripNumber[quadr].size()){
	StripNumber[quadr].push_back((*FirstHC)[j]->GetRadialID());
	StripEnergy[quadr].push_back(energy);
	if(data->VerboseLevel > 5){
	  cout<<"adding hit to new strip!"<<endl;
	}
      }      
    }//hit in active area
  }//firstHC entries
  if(SecondHC != NULL){
	 for(int j=0;j<SecondHC->entries();j++){
//      cout << (*FirstHC)[j]->GetSiSegmentID() << "	" << (*SecondHC)[j]->GetSiSegmentID() << endl;
		energy = (*SecondHC)[j]->GetEdep()/histoManager->GetEnergyUnit();
      int quadr = (*SecondHC)[j]->GetSiSegmentID();
      if(quadr >= 0 || data->WriteAll){ 
			if(quadr >= quads){
	  			cerr<<"error in hitcollection "<<i<<" hit "<<j<<" segmentid = "<<quadr<< " > " << quads <<endl;
	  			continue;
			}
		}
//	cout << "E: "<< energy << endl;
	Edet[quadr]+=energy;
    }
  }//there was a CDE Hit
  for(int j=0;j<quads;j++){
    if( (Etot[j]>data->DetectionEnergy) || (Edet[j]>data->DetectionEnergy) ){
      
		Annular cd;
      cd.ID(j);
      //cout << "sim strip " << StripNumber[j][0] << endl;
      cd.Time(data->EventTime);
      if( (RingNumber[j].size()==0) || (StripNumber[j].size()==0) ){
	if(data->VerboseLevel > 3)
	  cerr<<"Hit in CD sector "<<j<<" but no rings ("<<RingNumber[j].size()<<") or no strips ("<<StripNumber[j].size()<<") found, Edetector: "<<Edet[j] <<" deltaEdetector: "<<Etot[j] << endl;
	break;
      }
      if(RingNumber[j].size()==1)
	cd.SetRing(RingNumber[j][0], RingEnergy[j][0] + RandGauss::shoot(0, data->CDResolution/(2.*sqrt(2.*log(2.)))));
      else if(RingNumber[j].size()==2)
	cd.Set2Rings(RingNumber[j][0], RingEnergy[j][0] + RandGauss::shoot(0, data->CDResolution/(2.*sqrt(2.*log(2.)))),RingNumber[j][1], RingEnergy[j][1] + RandGauss::shoot(0, data->CDResolution/(2.*sqrt(2.*log(2.)))));
      else if(data->VerboseLevel > 3)
	cout << "Found more than 2 rings " <<RingNumber[j].size()<< endl;
      if(StripNumber[j].size()==1)
	cd.SetStrip(StripNumber[j][0], StripEnergy[j][0] + RandGauss::shoot(0, data->CDResolution/(2.*sqrt(2.*log(2.)))));
      else if(StripNumber[j].size()==2)
	cd.Set2Strips(StripNumber[j][0], StripEnergy[j][0] + RandGauss::shoot(0, data->CDResolution/(2.*sqrt(2.*log(2.)))),StripNumber[j][1], StripEnergy[j][1] + RandGauss::shoot(0, data->CDResolution/(2.*sqrt(2.*log(2.)))));
      else if(data->VerboseLevel > 3)
	cout << "Found more than 2 strips " <<StripNumber[j].size()<< endl;

//	cout << "Edet: " << Edet[j] << endl;

      if(Edet[j]>0)
	cd.SetEdet(Edet[j] + RandGauss::shoot(0, data->CDResolution/(2.*sqrt(2.*log(2.)))));

      if(forward)
	histoManager->FillFCDBranch(&cd);
      else
	histoManager->FillBCDBranch(&cd);
    }
  }
}  

//barrel
void EventAction::BD(bool forward, HitsCollection* FirstHC, HitsCollection* SecondHC){
  vector<int> StripNumber;
  vector<double> StripEnergy;
  StripEnergy.clear();
  //vector<double> StripPos;
  //StripPos.clear();
  double edet = 0.;
  double rear = 0.;
  int k;
  double energy;

  Int_t i_mod;

  i_mod = i - data->FirstBarrelHC; // i_mod runs now from 0 on


  if(data->ExperimentType < 2){
    cerr << "experiment type "<<data->ExperimentType<<" no barrel detector!" << endl;
    return;
  }
  int quadr = i_mod/data->BarrelDetectors;
  //cout << "starting barrel hit HC " << i <<" i_mod " << i_mod << " Quadrant " << quadr << endl;
  //cout << "entries " << FirstHC->entries() << endl;
  if(quadr >= data->BarrelDetectors){
    cerr<<"error in hitcollection "<<i<<" quadrant = "<<quadr<< " > " << data->BarrelDetectors <<endl;
    return;
  }
  for(int j=0;j<FirstHC->entries();j++){
    energy = (*FirstHC)[j]->GetEdep()/histoManager->GetEnergyUnit();
//  	cout << quadr << "	" << (*FirstHC)[j]->GetPos().phi()/3.1415*180.-90*(quadr==0)+90*(quadr==2)+180*((quadr==1)&&((*FirstHC)[j]->GetPos().phi()<0))-180*((quadr==1)&&((*FirstHC)[j]->GetPos().phi()>0)) << "	" << (*FirstHC)[j]->GetPosInDet().x() << endl;
//    	if (quadr==2) cout << (*FirstHC)[j]->GetPos().x() << "	" << (*FirstHC)[j]->GetPos().y() << "	" << (*FirstHC)[j]->GetPos().z() << endl;

	 if(energy==0.)
      continue;
    //cout << "energy " << energy << " strip " <<  (*FirstHC)[j]->GetAnnularID() << endl;
    if((*FirstHC)[j]->GetAnnularID() >= 0 || data->WriteAll){ //active area
      //check whether this strip was already hit
      //cout << "StripNumber.size() " << StripNumber.size() << endl;
      for(k=0; k<(int)StripNumber.size(); k++){
	if(StripNumber[k] == (*FirstHC)[j]->GetAnnularID()){
	  //cout << "same strip" << endl;
	  if(forward && data->settings->fBadStrips[0][quadr][(*FirstHC)[j]->GetAnnularID()]==0){
	    //cout << "same strip energy " << energy << " pos " << ((*FirstHC)[j]->GetPosInDet().x()+data->settings->BarrelActLength/2.)/data->settings->BarrelActLength << endl;
	    StripEnergy[k] += energy*((*FirstHC)[j]->GetPosInDet().x()+data->settings->BarrelActLength/2.)/data->settings->BarrelActLength + RandGauss::shoot(0, data->BarrelResolution/(2.*sqrt(2.*log(2.))));;
	  }
	  else if(!forward && data->settings->fBadStrips[1][quadr][(*FirstHC)[j]->GetAnnularID()]==0){
	    StripEnergy[k] += energy*(1-((*FirstHC)[j]->GetPosInDet().x()+data->settings->BarrelActLength/2.)/data->settings->BarrelActLength) + RandGauss::shoot(0, data->BarrelResolution/(2.*sqrt(2.*log(2.))));
	  }
	  break;
	}
      }
      //no matching channel was found
      if(k == (int) StripNumber.size()){
	//cout << "k "<< k << " strip " <<  (*FirstHC)[j]->GetAnnularID() << " energy " << energy << endl;
	if(forward && data->settings->fBadStrips[0][quadr][(*FirstHC)[j]->GetAnnularID()]==0){
	  //cout << "new strip energy " << energy << " pos " << ((*FirstHC)[j]->GetPosInDet().x()+data->settings->BarrelActLength/2.)/data->settings->BarrelActLength << endl;
	  StripNumber.push_back((*FirstHC)[j]->GetAnnularID());
	  StripEnergy.push_back(energy*((*FirstHC)[j]->GetPosInDet().x()+data->settings->BarrelActLength/2.)/data->settings->BarrelActLength + RandGauss::shoot(0, data->BarrelResolution/(2.*sqrt(2.*log(2.)))));
	  //StripPos.push_back((energy + RandGauss::shoot(0, data->BarrelResolution/(2.*sqrt(2.*log(2.)))))*((*FirstHC)[j]->GetPosInDet().x()+data->settings->BarrelActLength/2.)/data->settings->BarrelActLength/energy);
	}
	else if(!forward && data->settings->fBadStrips[1][quadr][(*FirstHC)[j]->GetAnnularID()]==0){
	  StripNumber.push_back((*FirstHC)[j]->GetAnnularID());
	  StripEnergy.push_back(energy*(1-((*FirstHC)[j]->GetPosInDet().x()+data->settings->BarrelActLength/2.)/data->settings->BarrelActLength) + RandGauss::shoot(0, data->BarrelResolution/(2.*sqrt(2.*log(2.)))));
	  //StripPos.push_back((energy + RandGauss::shoot(0, data->BarrelResolution/(2.*sqrt(2.*log(2.)))))*(1-((*FirstHC)[j]->GetPosInDet().x()+data->settings->BarrelActLength/2.)/data->settings->BarrelActLength)/energy);
	  //strip pos other way around for backward barrel
	}
      }
      rear+=energy;
            
    }//active area
  }

  if(SecondHC != NULL){
    for(int j=0;j<SecondHC->entries();j++){
      energy = (*SecondHC)[j]->GetEdep()/histoManager->GetEnergyUnit();
      if(quadr>=0 || data->WriteAll){ //active area
	edet+=energy;
      }
    }
  }
  if( (rear==0) && (edet==0) ){
    if( (data->VerboseLevel > 3) && (StripNumber.size()>0) )
      cerr <<"Error, found no hits in sector "<<quadr<<", but strips were hit (#strips = "<<StripNumber.size()<<", rear energy = "<<rear<<", Edet energy = "<<edet<<")"<<endl;
    return;
  }

  //smear energy resolution
  //cout << "sim rear " << rear << " edet " << edet << endl;
  rear += RandGauss::shoot(0, data->BarrelResolution/(2.*sqrt(2.*log(2.))));
  if(edet>0)
    edet += RandGauss::shoot(0, data->BarrelResolution/(2.*sqrt(2.*log(2.))));
  //cout << "det rear " << rear << " edet " << edet << endl;
  if(rear<0)
    rear=0.;
  if(edet<0)
    edet=0.;

  if( (rear > data->DetectionEnergy) || (edet > data->DetectionEnergy) ){  

    if(StripNumber.size()==0){
      if(data->VerboseLevel>2)
	cerr <<"Error, found no hits in sector "<<quadr<<"rear energy = "<<rear<<", Edet energy = "<<edet<<")"<<endl;
      return;
    }
    //for(int k=0; k<StripNumber.size();k++)
    // cout << "strip "<<k<<" " << StripNumber[k] << " quadrant " << quadr << endl;//testingcdposition
    Barrel bd;
    bd.ID(quadr);
    bd.Time(data->EventTime);
    bd.SetEdet(edet);
    bd.SetRear(rear);
    //if(!forward)
    //  cout << "backward pos in det " << StripPos[0] << " with energy " << rear << endl;
    
    bd.SetMult(StripNumber.size());
    if(StripNumber.size()==1){	
      //cout << "1 hit stripenergy " << StripEnergy[0] << endl;
      bd.SetStrip(StripNumber[0], StripEnergy[0]/rear);
      bd.SetMult(1);
    }
    else if(StripNumber.size()==2){
      bd.SetMult(2);
      //cout << "2 hits stripenergy " << StripEnergy[0] << " and " << StripEnergy[1] << endl;
      if(StripEnergy[0] > StripEnergy[1]){
	bd.Set2Strips(StripNumber[0], StripEnergy[0]/rear, StripNumber[1], StripEnergy[1]/rear,(abs(StripNumber[0]-StripNumber[1])));
      }
      else{
	bd.Set2Strips(StripNumber[1], StripEnergy[1]/rear, StripNumber[0], StripEnergy[0]/rear, (abs(StripNumber[0]-StripNumber[1])));   
      } 
    }
    else{
      if(data->VerboseLevel > 3)
	cout << "Found more than 2 strips " <<StripNumber.size()<< endl;
      bd.SetMult(3);
      int higheststrip = -1;
      int secondstrip = -1;
      int thirdstrip = -1;
      double highestenergy = 0.;
      double secondenergy = 0.;
      double thirdenergy = 0.;
      for(int k=0;k<StripNumber.size();k++){
	if( StripEnergy[k] > highestenergy ){
	  thirdstrip = secondstrip;
	  secondstrip = higheststrip;
	  higheststrip = StripNumber[k];
	  
	  thirdenergy = secondenergy;
	  secondenergy = highestenergy;
	  highestenergy = StripEnergy[k];
	}
	else if( StripEnergy[k] > secondenergy ){
	  thirdstrip = secondstrip;
	  secondstrip = StripNumber[k];
	  
	  thirdenergy = secondenergy;
	  secondenergy = StripEnergy[k];
	}
      }
      if(higheststrip>-1){

	bd.Set2Strips(higheststrip, highestenergy/rear, secondstrip, secondenergy/rear, (abs(higheststrip-secondstrip)==1));
      }
    }
    if(forward)
      histoManager->FillFBBranch(&bd);
    else
      histoManager->FillBBBranch(&bd);
  }

}

